# Neil J. Gunther's "Universal Scalability Model"
# M. Edward (Ed) Borasky
# licensed under GPL, just like PDQ.

# Usage:

# 1. Start up R. This has been tested under Linux and Windows and should work
#    other places as well.

# 2. Make sure you're in the right directory, then type

#    source("ScalabilityModel.R", echo=TRUE)

#    at the R "> " prompt.

# 3. Output is to PDF files. It will also write a log file, 
#    "ScalabilityModel.log".

# Define the function
UniversalScalabilityModel <- function(Users, Throughput, Output.PDF.File,
	trace=FALSE, ...) {

  # define function to be minimized
  to.be.optimized <- function(parameters) {
  
    # get and check parameters
    X1 <- parameters[1]
    if (X1 <= 0) return(1/0)
    sigma <- parameters[2]
    if (sigma < 0) return(1/0)
    if (sigma > 1) return(1/0)
    kappa <- parameters[3]
    if (kappa < 0) return(1/0)
  
    # now we have good parameters -- return sum of squared residuals
    fitted.Throughput <- X1*Users/(1 + sigma*(Users - 1) + 
      kappa*Users*(Users - 1))
    residuals <- Throughput - fitted.Throughput
    return(sum(abs(residuals)))
  }

  # now optimize
  opt <- optim(
    c(Throughput[1], 0, 0),
    to.be.optimized,
    method="Nelder-Mead"
  )
  print(opt) # show results
  if (opt$convergence != 0) return(w) # bail if it didn't converge
  X1 <- opt$par[1]
  sigma <- opt$par[2]
  kappa <- opt$par[3]
	
	# draw a smooth curve through the data
  n.curve.points <- 500
	z <- as.data.frame(seq(
		1, 
		max(Users), 
		by=max(Users)/n.curve.points
	))
  n.z.rows <- nrow(z)
	colnames(z) <- "Users"
  z$Throughput <- X1*z$Users/(1 + sigma*(z$Users - 1) + 
      kappa*z$Users*(z$Users - 1))
  z$Throughput[0] <- 0

	# compute plot labels
	Users.max <- sqrt(1/kappa - sigma/kappa)
  Throughput.max <- X1*Users.max/(1 + sigma*(Users.max - 1) + 
      kappa*Users.max*(Users.max - 1))
	Users.opt <- Throughput.max/(X1*(1 - sigma - kappa))

	# format subtitle with model parameters
	subtitle <- paste(
		"sigma = ", signif(sigma, digits=2),
		", kappa = ", signif(kappa, digits=2),
    ", X1 = ", signif(X1, digits=2),
	sep="")

	# format main title with results
	maintitle <- paste(
		"Universal Scalability Model: ", 
		Output.PDF.File,
		"\nUsers.max = ", signif(Users.max, digits=2),
		"\nThroughput.max = ", signif(Throughput.max, digits=2),
		", Users.opt = ", signif(Users.opt, digits=2),
	sep="")

	# make the picture
	pdf(Output.PDF.File, width=10, height=7.5) # open the plot file

	# plot the fitted curve
	plot(z$Users, z$Throughput, type="l", lwd=2, col="blue",
    ylim=c(0, Throughput.max),
		xlab="Users", 
		ylab="Throughput", 
		sub=subtitle,
		main=maintitle
	)
  grid(col="cornflowerblue")

	# plot the input capacity values as points
	points(Users, Throughput, col="darkgreen", pch=18)

  # and the bounds
  lines(c(1, Users.opt), 
    c(X1, Throughput.max), col="red", lwd=1)
  lines(c(1, z$Users[n.z.rows]), 
    c(Throughput.max, Throughput.max), col="red", lwd=1)

	# close plot file
	graphics.off()

	# print numeric results
	print(w)
	c <- try(confint(w), silent=TRUE)
	if (class(c)[1] != "try-error") {
		print(c)
	}

	# return the model
	return(opt)
}

# test cases from Gunther's spread sheet
#sink("ScalabilityModel.log") # divert output from console
z <- read.csv("cisco.csv")
w <- UniversalScalabilityModel(z$N, z$X, "cisco.pdf")
z <- read.csv("intel-was.csv")
w <- UniversalScalabilityModel(z$N, z$X, "intel-was.pdf")
z <- read.csv("nttest.csv")
w <- UniversalScalabilityModel(z$N, z$X, "nttest.pdf")
z <- read.csv("ntsql-4.csv")
w <- UniversalScalabilityModel(z$N, z$X, "ntsql-4.pdf")
z <- read.csv("ntsql-6.csv")
w <- UniversalScalabilityModel(z$N, z$X, "ntsql-6.pdf")
z <- read.csv("ntsql-8.csv")
w <- UniversalScalabilityModel(z$N, z$X, "ntsql-8.pdf")
z <- read.csv("sgi-numa.csv")
w <- UniversalScalabilityModel(z$N, z$X, "sgi-numa.pdf")
z <- read.csv("sun-e10k.csv")
w <- UniversalScalabilityModel(z$N, z$X, "sun-e10k.pdf")
z <- read.csv("spec-sdm.csv")
w <- UniversalScalabilityModel(z$N, z$X, "spec-sdm.pdf")
#sink() # return output to console
